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Abstract: We determine holographically 2-point correlators of gauge invariant operators 
with large conformal weights and entanglement entropy of strips for a time-dependent 
anisotropic 5-dimensional asymptotically anti-de Sitter spacetime. At the early stage of 
evolution where geodesics and extremal surfaces can extend beyond the apparent horizon 
all observables vary substantially from their thermal value, but thermalize rapidly. At late 
times we recover quasi-normal ringing of correlators and holographic entanglement entropy 
around their thermal values, as expected on general grounds. We check the behaviour of 
holographic entanglement entropy and correlators as function of the separation length of 
the strip and find agreement with the exact expressions derived in the small and large 
temperature limits. 
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1 Introduction 


Far from equilibrium dynamics and thermalization of strongly coupled systems has attracted 
a lot of attention in the past decade. The reason for this is twofold. On the one hand 
experiments at RHIC and the LHC revealed that the quark gluon plasma created in heavy 
ion collisions behaves as a strongly coupled liquid that thermalizes extremely fast [1-4]. 
On the other hand, in condensed matter experiments it is now possible to drive an isolated 
system to a far from equilibrium state by a quantum quench, i.e. a control parameter of 
the system is varied rapidly [5]. 

One powerful tool to study strongly coupled systems out of equilibrium is the gauge 
gravity duality where classical supergravity in d + 1 dimensions is dual to a d-dimensional 
strongly coupled large-held theory that loosely speaking lives on the boundary of the 
gravity theory [6-8]. The thermalization process of held theories with a conformal held 
theory hxed point in the ultraviolet is then mapped to black hole (or black brane) formation 
in asymptotically Anti-de Sitter (AdS) space [9]. 

Particularly useful theoretical observables to monitor the thermalization process are 
the vacuum expectation value of the stress-energy tensor, correlation functions of gauge 
invariant operators and entanglement entropy (EE). The former two are easy to dehne on 
the held theory side and, given some standard assumptions that we recall below, easy to 
calculate on the gravity side. We recall now relevant aspects of EE, see [10] for more details 
and references. 

On the held theory side EE is dehned in the following way. Dividing a system into 
two subsystems A and B the EE Sa of the subsystem A is dehned as the von Neumann 
entropy of the reduced density matrix obtained by tracing out the degrees of freedom of 
the subsystem B. 

On the gravity side, the holographic entanglement entropy (HEE) is dehned as the area 
of a minimal surface extending from some predehned surface A on the boundary into the 
bulk [11, 12]. For time-dependent backgrounds the minimal surface has to be replaced by 
an extremal surface [13]. It is noteworthy — both for conceptual and for technical reasons 
— that the relevant surfaces used to determine HEE can extend beyond the apparent 
horizon. Conceptually, this implies that HEE provides a ‘classical’ way (on the gravity 
side) to extract (quantum) information from the region beyond the horizon. Technically, 
this requires the numerical determination of (part of) the spacetime region beyond the 
apparent horizon, which can be a challenge since that region ends in a singularity. 

EE is notoriously hard to compute in quantum held theories; so far this is only possible 
in highly symmetric theories such as (relativistic) conformal held theories [14-16] or Galilean 
conformal held theories [17]. In a seminal work, Calabrese and Cardy [18] were able to 
compute the time evolution of the EE after a quench in a two dimensional conformal held 
theory and in the Ising spin chain model in a transverse magnetic held. In both cases they 
hnd that for an entangling interval of length I, the EE increases linearly with time until 
t ^ lj2 after which it saturates. The linear scaling with time and the crossover at t ~ 1/2 
can be understood in terms of entangled quasiparticles pairs emitted from the initial state 
and is therefore expected to hold for a wider class of systems. 
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The holographic duality offers a playground where one can study various different 
systems and search for universal and novel properties of HEE in dynamical situations. 

The simplest example where one can study the analog of quenches in the holographic 
setup are spacetimes where thin shells collapse to form a black hole, first utilised in [19]. 
The behaviour of the EE in these setups has been studied extensively [19-25] and indeed 
shows universal behaviour consistent with the findings of Cardy and Calabrese. Namely 
the initial short early time epoch is followed by a long period where the EE grows linearly 
with time and is independent of the entangling boundary region. This was first worked out 
in the Vaidya spacetime [19, 25, 26] where the shell is composed out of null dust and later 
generalized to matter with arbitrary equation of states [27, 28]. The linear scaling is even 
present in geometries with Lifshitz scaling and hyper scaling violation [29, 30]. In addition, 
in the above works the HEE is a monotonically increasing function that saturates to its 
equilibrium value from below. Note that these quenches are thermal quenches because the 
end state is a thermal state. 

Another way of studying thermal quenches is by turning on a radially collapsing scalar 
field that forms a black hole. This situation is more complicated because in order to obtain 
the geometry from which one can then obtain the HEE, Einsteins equations have to be 
solved numerically [31-33]. In [32, 34] this was done for a radial collapsing massless scalar 
field in global AdS where the scalar field can have many bounces between the boundary 
and the center of AdS before a black hole forms resulting in a periodic behaviour of the 
HEE. 

In [33] a massive scalar field dual to a massive fermionic operator was turned on, 
treating the quench as a perturbation on the static spacetime. In this setup the HEE is 
also not monotonic and in some cases approaches the equilibrium value from above. This 
reveals a qualitative difference of EE to the thermal entropy which, on general grounds, 
must be monotonically growing in a closed system. 

An additional motivation to study HEE comes from the question how to measure 
entropy production in (holographic models for) heavy ion collisions [35]. Within the 
gauge/gravity duality the entropy of the (stationary) black hole corresponds to the en¬ 
tropy of the field theory. However, in time-dependent backgrounds entropy as defined from 
the area of the apparent horizon is ambiguous because it depends on the choice of time 
slicing. By contrast the definition of the HEE is unique and therefore may serve as an 
alternative measure for entropy production. 

So far nearly all studies of HEE have relied on the simplifying assumption of a spheri¬ 
cally symmetric spacetime (see [36] for a notable exception). In the paper at hand we drop 
this assumption and investigate the effect of anisotropies on the evolution of the HEE. 

The paper is organized at follows. In section 2 we introduce the anisotropic spacetime 
we will use, recall general aspects of 2-point correlators and simplify the determination of 
HEE to a geodesic problem in an auxiliary spacetime. In section 3 we discuss the numerical 
strategy that we used, relegating details to the appendices. In section 4 we present results 
for the background geometry, holographic stress tensor, 2-point correlators and HEE. In 
section 5 we analyze the late time behaviour of correlators and HEE and relate it to the 
lowest lying quasinormal mode. In section 6 we conclude with a brief summary of our 
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results and some possible next steps. In appendix A we comment on the spectral method 
and other numerical routines used to solve the 5-dimensional vacuum Einstein equations. 
In appendix B we explain the relaxation code used to determine geodesics and extremal 
surfaces. 

Before starting we mention some of our conventions. We use mostly -t- signature and 
set the AdS radius and final black brane mass to unity. 

2 Theoretical setup 

2.1 Anisotropic asymptotically AdSs spacetimes 

In this section we review the most important details of the model first introduced in [37] 
and studied further in [38-41]. 

The 5-dimensional bulk metric that introduces anisotropy between the longitudinal 
and transverse directions with an 0(2) rotational invariance in the transverse plane can be 
written conveniently in Eddington-Finkelstein coordinates, 


= —A{r,v)dv^ + 2<lrdv + T,‘^{r,v){e -|-dx^^ (2.1) 

where the functions A, B and S only depend on the holographic coordinate r and (ad¬ 
vanced) time V. In this coordinate system the vacuum Einstein equations read 

0 = S(S)'-h 2S'i] - 2^2 (2.2a) 

0 = 2'L{B)' + 3{T,'B + B't) (2.2b) 

d = A"+ 3B'B-l2T.'t/T?+ d (2.2c) 

0 = 2S - A'S-h (2.2d) 

0 = 2S" -h {B'fT. (2.2e) 

where prime denotes radial derivative and dot time derivative, viz. 

h' = drh h = dyh + ^Adrh (2.3) 


for any function h{r,v). 

The Einstein equations (2.2) have to be solved for special initial conditions and appro¬ 
priate boundary conditions. There are, at least, two ways to create a far from equilibrium 
state. On the one hand one can turn on a time dependent anisotropy function at the bound¬ 
ary B{r = oo,v) = BQ{t) as in the original works of [37, 38] and let the system evolve. In 
this case the boundary metric is curved and the conformal anomaly is present [42]. On the 
other hand one can specify the initial state in the absence of external sources by specifying 
the metric in the bulk on the initial time slice [39] with a flat boundary geometry. For 
simplicity, in the following we will study the setup where the boundary metric is flat and 
time independent. 
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Requiring that the spacetime is asymptotically AdSs, Einstein’s equations in the near 
boundary expansion (r —?• oo) are solved by^ 






B = hW + ^ + 0(r-«) 

p^ pO '■ '' 


(2.4a) 

(2.4b) 

(2.4c) 


The coefScients in the asymptotic expansion determine the expectation value of the stress 
energy tensor in the dual field theory [43] 




iV2 

^ ^ c 

2tt'^ 


diag[T, P||(t), 


P±it), P±{t)] 


(2.5) 


where 





1 

4 


04 — 


P±{t) 


--04 + hi{t). 


( 2 . 6 ) 


In order to determine the fourth order coefficients one needs to solve Einstein’s equations 
numerically for some initial conditions. When doing the actual calculation we work with 
the inverse radial coordinate z = 1/r so that the boundary is located at z = 0 . 

For our initial data we follow [40, 41] and choose for the anisotropy function on the 
initial time slice 


B{r,vo) = 4 exp 



(2.7) 


with f3 = 6.6, ro = 4 and w = 1. In addition the initial conditions have to be supplemented 
with a value for the coefficient 04 which sets the energy density of the initial state, for which 
we take 04 = —1, corresponding to an equilibrium temperature T = I/tt, as we now recall. 

At late times we expect isotropization, B = 0. In that case we recover the usual static 
AdS black brane solution as follows. Solving (2.2e) and using residual gauge transformations 
yields S = r. This implies S' = 1 and S = | A. Solving (2.2a) then yields A = (1 —1/r^), 

where we fixed the integration constant such that 04 = —1. [The other equations are either 
trivial, (2.2b) and (2.2d), or redundant, (2.2c).] The result for A is the usual Killing norm 
for the static AdS black brane. Surface gravity is given by k = ^ 1 = 2 so that the 

Hawking temperature is T = K/{2Tr) = I/tt. 

In the generic anisotropic case, R 7 ^ 0, we solve the Einstein equations (2.2) numerically 
for the initial conditions (2.7). In this background we then study the evolution of 2-point 
correlation functions for operators of large conformal weights and the HEE. This in turn 
requires us to determine the background sufficiently far beyond the apparent horizon. In 
section 3.1 below we discuss the numerical implementation. 


^We fix a residual gauge freedom that preserves the form (2.1), r —^ r + f{v), by demanding that the 
first subleading term in the function A falls off like 0{r~^). Note that in several numerical simulations in 
the literature the same freedom is fixed by placing the apparent horizon at a specific value of the radial 
coordinate r. 


- 5 - 









2.2 2-point correlators 


The equal time 2-point function for an operator of large conformal weight A can be com¬ 
puted via a path integral as [44, 45] 

{0{t, x)0{t, /)) = J ~ ~ (2.8) 

geodesics 

where the integral is a sum over all possible paths with endpoints at {t, x') and (t, x) 
and C{V) is the proper length of the path. The hrst approximation neglects perturbative 
corrections and is the so called geodesic approximation, which holds in the limit when the 
conformal weight of the operator is large. The conformal weight effectively plays the role 
of 1/h in usual perturbative expansions of path integrals. Then it can be shown that the 
sum over all paths reduces to a sum over all geodesics where Lg denotes the length of the 
corresponding geodesic. To leading order only the geodesic with the smallest value of Lg 
contributes, whose length we denote by L, which explains the second approximation^. It 
neglects instanton corrections. 

However, the length of the geodesic has a divergence originating from the asymptotically 
AdS boundary and therefore needs to be renormalized. We choose to subtract the length 
of a geodesic in the static black brane background, which we denote by Ttherm- In terms of 
the renormalized length 5L = L — Ttherm the 2-point function becomes 

{0{t,x)0{t,j^)) r. . (2.9) 

This means that we can obtain the time evolution of 2-point functions by looking at spacelike 
geodesics that are anchored at the boundary at hxed separation I and calculating their 
length at different times. Due to the anisotropy in the system we only solve for the subset 
of correlation functions that are either separated in the longitudinal direction or in the 
transverse directions. 

To this end we let all the coordinates depend on one parameter a, which lies in the 
interval a G [—cr^, CTm]- To obtain the lengths of the geodesics we have to solve the geodesic 
equation for the two subspaces given by the line elements 

dis\ = —Adv^—dz dv + T? d.x\ ( 2 . 10 ) 

dsf, = -Adv"^ - ^dzdv + Y^e-"^^ dxl. ( 2 . 11 ) 

II II 

For the separation in the transverse direction the geodesics end at (u(±fTm) = t, 
x_Lj(±(Tm) = diXol2, x_l 2 (±cTm) = 0, X|| (icTm) = 0), where t is the boundary time. Simi¬ 
larly, for the longitudinal separation we take (u(±cJm) = t, iC||(ic’'m) = ixo/2, x±{dzam) = 
0). With this choice of boundary conditions the lengths of the geodesics in the background 

^ For a comparison of the 2-point correlation function obtained by using the “extrapolate” dictionary 
and the geodesic approximation in AdSa Vaidya spacetime see [46]. 
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( 2 . 1 ) are given by 


/ ^m, / O 

daJ—A{v')‘^ -(2.12) 

/ Cm I 2 

day —A{v')^ - ^z'v' + (2.13) 

where prime denotes the derivative with respect to a. 

It is important to point out that we can only study geodesics after some advanced 
time V > Umin with boundary separations below a maximal separation I < /^ax- This comes 
from the fact that by solving Einstein’s equations numerically we have to choose a finite 
computational domain. Also, by specifying the initial state in the entire bulk on the initial 
time slice the advanced time interval at our disposal is G [r;o,oo]. As the geodesics 
reach into the bulk they bend back in advanced time leaving the computational domain for 
advanced times v < Umin as well as extending too far into the bulk for separations I > /max- 
In section 3.2 below we discuss how to solve the geodesic equations numerically. 

2.3 Holographic entanglement entropy 

In time dependent systems the covariant HEE [13] for some boundary region A is obtained 
by extremizing the 3-surface functional 

.4 = /dV^det(|^^9,„) (2.14) 

that ends on the boundary surface A. In the dual field theory the EE is then conjectured 
to be given by [11, 13, 47] 

5be = ^. (2.15) 

Usually the boundary regions of interest are either a sphere or a strip that has finite 
extent in one direction and infinite extent in the other two directions. In spacetimes with 
spherical symmetry in the three spatial dimensions the problem of finding the extremal 
area functional (2.14) effectively reduces to finding geodesics. In our case where spherical 
symmetry is broken this is not the case anymore. For example, finding the extremal area 
for a spherical boundary region would require to solve nonlinear coupled partial differential 
equations. However, in the case of a strip with finite extent either in the transverse or 
longitudinal direction it is possible to reduce the problem to finding geodesics in a suitable 
auxiliary spacetime, as we now demonstrate. 

We introduce two scalar fields (/>j(x“) and write the line element as 

ds^ = dx^ dx'^ = hag dx“ dx^ + (fl dx 2 + 4>2 (2.16) 

where hag is a 3-dimensional metric with coordinates {v,r,xi) where xi represents the 
coordinate we choose to have finite spatial extent, i.e. either xy or one of x_l. The remaining 
(non-compact) coordinates are then denoted by X 2 , X 3 , which we choose to be two of 
our three world-volume coordinates; the third one is denoted by a. Parametrizing the 


- 7 - 







3-dimensional coordinates as x" = {v{a),r{a),xi{a)), the area functional (2.14) can be 
written as 

A = jdx^j dx 2 j do\j ■ (2.17) 

Performing the integration over the Killing coordinates X 2 and X 3 yields a (possibly infinite) 
constant volume factor through which we are going to divide. Thus, instead of calculating 
HEE we calculate a HEE density per Killing volume. The problem of extremizing the 
3-surface corresponding to a boundary region A of strip-topology is then reduced to a 
1 -dimensional problem. 

In fact, from the expression (2.17) on can see that the problem of finding the extremal 
3-surfaces reduces to finding geodesics of the conformal metric 

ds = hais dx°‘ dx^ = (t){(t>2hai3 dx" dx^ . (2.18) 

The 3-dimensional conformal metrics for separation in the transverse and longitudinal di¬ 
rections for which we have to solve the geodesic equation in our case are given by 

ds5_ = —Adv^+ 2drdv dx\^ (2.19a) 

dsj = Adv"^+ 2 dr dv + Y^e-"^^ dxf) . (2.19b) 

3 Numerical implementation 

3.1 Einstein equations 

We solve the Einstein equations (2.2) using pseudo spectral methods as described in detail 
in [41], with the only difference that we do not fix the location of the apparent horizon, 
see appendix A for some details. Not fixing the apparent horizon facilitates the study 
of geodesics and extremal surfaces that reach behind the apparent horizon, which is of 
relevance for 2-point functions and HEE. For that reason we want a large computational 
domain in the holographic coordinate z = 1/r. In all the computations we took z G [0, 1.6] 
with the final position of the horizon located at z = 1. For the time evolution it is sufficient 
to use a fourth order Runge Kutta method with time steps 5t = 10“^. All the computations 
were done with the open source software GNU Octave [48]. 

3.2 Geodesics 

To compute 2-point functions we need to find curves of extremal length in a curved space- 
time whose endpoints reside on fixed positions on the boundary of that spacetime. These 
curves are solutions to the geodesic equation subject to boundary conditions at the end¬ 
points. For numerical reasons it turns out to be convenient to use a non-affine parameter 
a , where r = r(cT) is the usual affine parametrization, = 1 - In terms of a the 

geodesic equation reads 

= -JX^ (3.1) 

where X^ = and J = denotes the Jacobian which originates from the change 

in parametrization. This form of the geodesic equation gives us the freedom to choose 
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parametrizations resulting in better convergence behaviour of the relaxation algorithm than 
the affine parametrization does. In physical terms the right hand side in (3.1) introduces a 
fictitious viscous force that enhances numerical convergence. 

In our case the geodesic equation is given by a set of three coupled nonlinear ODEs 
of second order for the geodesic coordinates V,Z,X. This set of equations can be reduced 
to a set of six first order equations in terms of the geodesic coordinates and their first 
derivatives: 


pv = V (3.2a) 

pz = Z (3.2b) 

px = X (3.2c) 

Pv+ ^^vvPv+ xxPx = —Jpv (3.2d) 

Pz + Pv + PvPZ + Pz + Px = -Jpz (3.2e) 

Px + Zr^vxPvPx + zxPzPx = -Jpx (3.2f) 


This set of equations is a two-point boundary value problem, which is usually either 
solved with shooting methods or relaxation methods [49]. We do not shoot but relax. The 
geodesics obtained in appendix B serve as the initial guess for the relaxation algorithm. 
Since we are interested in one-parameter families of geodesics (evaluated at different con¬ 
stant time slices) we can take the solution for the family member as initial guess for 
the (n -|- 1)®* family member. More details of our implementation of the relaxation method 
are described in appendix B. 

3.3 Extremal surfaces 

The same method as above works also for HEE. Namely, for our problem at hand the 
evaluation of extremal surfaces is reduced to the evaluation of geodesics in some auxiliary 
spacetime, as we have shown in section 2.3. 

4 Results 

In this section we display and discuss our main results. In all figures where a time-axis 
is plotted we measure the boundary time t or the bulk advanced time v in units of the 
temperature of the final black brane, T = l/vr. The separation length I of 2-point functions 
and HEE and the corresponding boundary coordinates are given in units of T as well. To 
make the approach to thermal equilibrium most transparent we use normalized quantities 
for the geodesic length L^en and HEE S^en defined by 

Lren=^^^^ (4.1a) 

-^th 

Sren = (4.1b) 

where L (S') is the unrenormalized length (HEE) and Lth (Sth) is the corresponding thermal 
value. 
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Figure I. Left: anisotropy function B{r,v). Right: transverse and longitudinal pressure. 


4.1 Background geometry and holographic stress tensor 

Figure I displays the most salient features of the background geometry. The left figure plots 
the anisotropy function B(r,v) and displays the regions outside and inside the apparent 
horizon, as well as the event horizon. The black lines depict a null congruence of geodesics 
close to the event horizon to exhibit their ingoing/outgoing nature. The right figure plots 
transversal and longitudinal pressures as function of boundary time t. Note the quick 
thermalization of the pressure components. 

4.2 2-point correlators 

As we have noted before geodesics can extend beyond the apparent horizon. This is made 
explicit in Fig. II where the blue curve serves as our initial guess for the relaxation code. 
The red geodesics at late times approach the apparent horizon without crossing it. At 
sufficiently early times (and sufficiently large separation) the geodesics cross the apparent 
horizon, an example of which is depicted by the cyan curve. 

The evolution of the renormalized lengths in the transverse and longitudinal direc¬ 
tions for different separations are depicted in Fig. III. Depending on the separation the 
2-point functions start at t = tmin which is the time when the geodesics extend beyond the 
computational domain. 

The first observation is that the transverse and longitudinal directions oscillate out 
of phase as shown in Fig. IV. The same feature is seen in the transverse and longitudi¬ 
nal pressure. By comparing the thermalization times of the one point functions, i.e. the 
expectation value of the stress energy tensor with the 2-point functions we see that the 
2-point functions thermalize later as expected. Also, the thermalization time increases if 
the boundary separation is increased. 
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Figure II. The green (black) line indicates the 2-position of the apparent (event) horizon; the 
dark blue curve is the Poincare patch AdS geodesic we use to initialize the simulation; red curves 
are geodesics with different boundary separation probing the thermal regime (none of them crosses 
the apparent horizon); the cyan curve in the left part of each plot is a geodesic which probes the 
non-thermal regime and reaches beyond event and apparent horizon. Isometric view (left) and view 
in a;-direction (right). 


Longitudinal Separation Transverse Separation 




tT tT 


Figure III. Renormalized length of geodesics for different separations in longitudinal and transverse 
directions. 


4.3 Holographic entanglement entropy 

The extremal surface equations — which we mapped to geodesic equations in an auxiliary 
spacetime — are solved again by a relaxation method. We observe the same qualitative 
features as for geodesics in Fig. II above: at early times extremal surfaces can extend 
beyond the apparent horizon, while at sufficiently late times they approach it from the 
outside without crossing. However, there are also notable differences to geodesics, which 
we discuss now. 
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Longitudinal vs. Transverse Separation 



tT 


Figure IV. Comparison of longitudinal and transverse geodesic lengths for the same boundary 
separation. 


Longitudinal Separation Transverse Separation 




tT tT 


Figure V. Longitudinal and transverse HEE for different separations. 


As can be seen from Fig. IX in appendix B.2 conformal geodesics reach much further 
into the bulk compared to the pure AdS case. Therefore the boundary separations we 
can study for the HEE are smaller compared to the 2-point functions. This is also the 
reason why for the same boundary separation the HEE reaches equilibrium later as the 
2-point functions. For the same boundary separation and at the same boundary time 
conformal geodesics reach deeper into the bulk and further back in time and therefore are 
more sensitive to out of equilibrium effects which are most pronounced at early times. In 
addition the shape of the curves differ from the 2-point functions with the oscillations less 
pronounced. We exhibit these features now in some plots. 

Figure V plots HEE for different separations in logitudinal and transverse directions. 
Comparison with Fig. HI shows that the oscillations are less pronounced for HEE. 

Figure VI plots HEE for a fixed separation in longitudinal and transverse directions. 
Again the behaviour of the curves is out of phase, in the sense that maxima of one curve cor¬ 
respond to minima of the other. Comparison with Fig. IV shows again that the oscillations 
are less pronounced for HEE. 
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Longitudinal vs. Transverse Separation 



tT 


Figure VI. Longitudinal and transverse HEE for same separation. 


5 Late time behaviour and quasinormal modes 

After the early far from equilibrium phase the geometry relaxes to the static Schwarzschild 
black brane solution. As noted in [40, 41] the anisotropy of the system is exponentially 
damped and at sufficiently late times one enters the linearized regime. In this regime 
the approach to equilibrium is accurately described by the lowest lying quasinormal mode 
(QNM) which characterizes the response of the system to infinitesimal metric perturbation. 
In the case at hand the relevant channel for the gravitational fluctuations is the spin two 
symmetry channel which coincides with the fluctuations of a massless scalar field in the 
static black brane geometry. The asymptotic response of the pressure anisotropy then 
takes the form 


b^it) ~ Re (5.1) 

with the lowest QNM given by [50, 51] 

^ = ±3.119452 - 2.746676 T (5.2) 

On the field theory side QNMs appear as poles in the retarded Green function [51-54]. It 
is therefore expected that also the late time behaviour of the correlation functions obtained 
in the previous section is described by the lowest QNM. We now show that this is indeed 
the case. 

In figure VII (left) we plot the renormalized geodesic length multiplied with the imag¬ 
inary part of the lowest QNM for transverse and longitudinal separations. 

One clearly sees that after a short period of time the evolution of the correlator is accu¬ 
rately described by the ringdown of the black brane with constant amplitude and frequency. 
The connection between the late time behaviour of correlation functions and QNMs was 
previously also observed in [55-57]. 

It turns out that HEE also follows this pattern. In [58] a connection between QNMs 
and the behaviour of the HEE was found. From linearised Einstein equations one can derive 
a differential equation for the first order correction /S.Sa of the HEE describing its change 
when a given ground state is excited. By imposing infalling boundary conditions at the 
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Figure VII. Left: Renormalized geodesic length for longitudinal (red) and transverse (blue) sep¬ 
aration for IT = 0.32 multiplied by the imaginary part of the lowest QNM. Right: Renormalized 
HEE for the same parameters as on the left. 




Figure VIII. Left: Renormalized length of geodesics as a function of their boundary separation in 
transverse directions for different fixed boundary times t = 0.5, 1, 1.5 (endpoints from left to right). 
The curves terminate when the geodesics leave the computational domain. The black dashed line 
shows the thermal limit. Right: HEE (blue) and geodesic length (red) in the thermal (solid) and 
zero temperature (dashed) limit. 


horizon one obtains a QNM dispersion relation putting a constraint on HEE. With our 
numerical solution we can demonstrate that the late time behaviour of HEE indeed follows 
the QNM ringdown even without imposing infalling boundary conditions. In Fig. VH (right) 
we show the HEE multiplied with for the infinite strip with finite separation 

in longitudinal and transverse direction. As for the correlation function, at late times, the 
HEE shows quasinormal ringing with constant amplitude and frequency. These oscillations 
show that HEE must not approach its thermal value from below but rather shows oscillatory 
behaviour around its thermal value. 

To conclude this section we finally study the departure of the length of the geodesics 
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and HEE from equilibrium for different times as a function of the boundary separation. 
This time we normalized the length of geodesics by subtracting a cutoff dependent piece. 
The case for the 2-point function with separation in the transverse direction is displayed in 
Fig. VIII. Out of equilibrium effects manifest themselves as oscillations around the thermal 
value. The curves terminate when the geodesics leave the computational domain, so for 
early times we only have access to rather small boundary separations. The same effect is 
seen for HEE. 

In the thermal limit the scaling for the geodesic length at small and large boundary 
separation is dictated by conformal symmetry and is proportional to 21og(//2) and 21 
respectively. At large separation the HEE also scales linearly with the separation length, 
whereas at small separation it is proportional to All these results agree precisely with 
the perturbative expressions derived in the limits of small and large temperatures [47, 59]. 
Our numerical results are shown in Fig. VHI where we also plot the corresponding zero 
temperature results, which coincide with the thermal curves for small separations. 

6 Conclusions 

In the paper at hand we have studied the time evolution of 2-point correlation functions in 
the geodesic approximation and holographic entanglement entropy in an anisotropic AA = 4 
SYM plasma. 

To obtain the 2-point correlation function for separation in transverse and longitudinal 
direction we solved for the length of the geodesics in the anisotropic background (2.1) by 
using a relaxation method with the zero temperature geodesic as an initial guess. At the 
n-th step the solution of the (n — l)-st step is used as initial guess. 

Choosing an infinite strip with finite separation either in the transverse or longitudinal 
direction and using the symmetries of the system the problem of finding minimal surfaces 
was reduced to finding geodesics in an auxiliary conformally related spacetime and the same 
strategy as for the 2-point correlation function goes through. 

The 2-point correlation functions as well as holographic entanglement entropy show os¬ 
cillatory behaviour around their thermal value, where transverse and longitudinal directions 
oscillate out of phase. Both quantities approach their equilibrium value by exponentially 
damped oscillations and after the initial far from equilibrium epoch the thermalization pro¬ 
cess is accurately described by the lowest quasinormal mode. That the late time behaviour 
of holographic entanglement entropy is captured by the ring down of the black brane geome¬ 
try is one of the main results of this work and to our knowledge the first explicit verification 
of the connection between QNMs and holographic entanglement entropy found in [58]. 

The methods developed and applied in this work can be generalized to other time- 
dependent asymptotically anti-de Sitter backgrounds of interest, such as colliding shockwave 
backgrounds [60, 61], which have been constructed numerically in the past few years [62- 
64]. Another interesting generalization would be the consideration of compact entangling 
regions, e.g. with spherical topology to check the (in-)dependence of our results on the form 
of the entangling region. 
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A Spectral method 

in this appendix we give some details on how we numerically solve the Einstein equations 
(2.2). We work with the inverse radial coordinate z = 1/r such that the boundary is 
located at z = 0. Due to asymptotic AdS boundary conditions, the metric functions A 
and S diverge as z —>■ 0. It is numerically favourable to dehne new functions with the 
known divergent pieces removed and rescale them with appropriate powers of z so that the 
resulting functions are hnite or vanish as z —>■ 0 ; the precise boundary conditions on the 
new functions are presented in (A. 6 ) below. This leads us to the following held redehnitions 


A{z,v) ^ + zA{z,v) B{z,v) z^B{z,v) 

^{z,v) ^+ z2S(z,n) S(z,n) ^ + ^z^S(z, n) 

S(z, v) ^ ^ + 2 S(z, v) B{z, v) -2z^B{z, v). 


(A. la) 
(A.lb) 
(A.lc) 


The redehned anisotropy function B allows to extract 64 (t) simply from the boundary value 


of B' 


64 (t) = B' (z = 0,v = t) 


(A.2) 


where B' = dzB. In terms of the redehned helds the hrst four Einstein equations (2.2) can 
be rewritten in the form 



(A.3b) 


(A.3a) 


(A.3c) 


A" + ^A' + 4^ = JA 
z z^ 


(A.3d) 
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with the source functions given by 


JS = 

--zB^ -3z‘^BB' - -z^B'^ 

2 2 

(A.4a) 

= 

2(5E + 2z3s2 + zE') 

-|- z^T, 

(A.4b) 

3b ~ 

3{l + z^t){3B + zB') 

(A.4c) 

8 z 2 (l + ^35]) 

6(E(4 + z^T.) + z^B{l + z^T.f{3B + zB’) + zE' - zS(l - 2^^^ - z^s')) 

JA = 

(A.4d) 

(z + z4E)2 


The relation between dot-derivative and time-derivative, originally given in (2.3), turns into 


B = + \b' + ^B + ^z^AB + ]z^B'A. (A.5) 

2 A Az A A 

The boundary conditions for the redefined fields read 

T,{z = 0,n) = 0 T;'{z = 0,n) = 0 (A.6a) 

'E{z = 0,n) = 04 (A.6b) 

B{z = 0,n) = B'{z = 0,n) (A.6c) 

A(z = 0,v) = 0 A'{z = 0, n) = 04 (A.6d) 


where we set 04 = —1. On the initial time slice the boundary condition (A. 6 c) is computed 
from the initial conditions 


B{z, vq) = Pzexp 


{z - zof- 


(A.7) 


where we set j3 = 6 . 6 , zq = 1/4 and a; = 1. 

For a given set of initial and boundary data the system of equations (A.3) allows for 
the following solution strategy: 


1. With the initial conditions (A.7) and the boundary conditions (A. 6 a) we first solve 
(A.3a) for S on the initial time slice. 

2. For given E and the boundary condition (A. 6 b) we solve (A.3b) for E. 

3. Having B, E and E we next solve (A.3c) for B using the boundary condition (A. 6 c). 

4. With B, E, E and B and the boundary conditions in (A. 6 d) we solve (A.3d) for A. 

5. Finally we integrate (A.5) to get the new B on the next time slice and repeat the 
whole procedure all over again. 


For constant v we solve each of the equations (A.3) with a pseudo-spectral method [65] 
where the A^ -|- 1 grid points Ui on an intervall [a, b] are located at 


= 2 ((“ + b) + (a-b) cos{iTT/N)) 


i = l,...,N + 1. 


(A. 8 ) 


17 - 








In our computations we choose Zi G [0,1.6] and = 60. These points can be used to 
construct the entries of the spectral differentiation matrix Dij [66]: 


Doo = 

2 V2 + 1 

DjsfN = 

2N^ + 1 

6 

6 

— 


J = 1, ■ 

..,iV-l 

2(1 - z]) 

ci 

Dij = 

j / J 

. 

II 

Cj {Zi - Zj) 


(A.9) 
(A. 10) 

(A.ll) 


where cq = cjv = 2 and otherwise c* = 1. The derivative of a function f{zi) = fi on the 
spectral grid points is obtained by multiplication with this diffentiation matrix 


n = . 


(A.12) 


This allows to turn each of the equations (A.3) into a system of linear equations. For 
example the second equation in (A.3) translates to 




where the matrix Lij is given by 


Lij = Dij + diag 


2z‘^{3T, + zS') 


and the source vector by 


{jt)i — 


1 + z^T, j i] 

2(5S + 2z3s 2 + zS') 


(A. 13) 


(A. 14) 


(A.15) 


The boundary condition Si = 04 is implemented by setting (j^)! = 04 and Lij = 5ij. The 
solution vector Sj is then obtained by multiplying the inverse of Lij with the source vector 


Si — [Lij] ij±)j ■ 


(A. 16) 


We do this with the OCTAVE routine linsolve. The equations for Sj, Ai and Bi can be 
solved in the same way. 

To advance the solution for B to the next time slice it is sufficient to use a simple 
fourth order Runge Kutta method [49] 

B{z, V + 5v) = B{z, v) + 6v (^^ki + ^k 2 + 5^3 + ^kij. (A.17) 

where the koefficients ki are given by 


ki = dj;B{z,v) (A.18) 

k 2 = dv{B{z,v) + ^ki) (A.19) 

ks = dy{B{z,v) + 5 /^ 2 ) (A.20) 

ki = dv{B{z, v) + ks) (A.21) 


and dyB is computed from (A.5). In our simulations we use a stepsize of 6v ~ 0.01 which 
is sufficient to achieve a stable time evolution. 
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B Relaxation method 


In relaxation methods differential equations are replaced by approximate finite difference 
equations (FDEs) on a discrete set of points. The solution is determined by starting with 
an inital guess and improving it iteratively. In this iterative procedure the result is said to 
relax to the true solution. 

First we define a grid ai = hi of equidistant spacing h = with i = ... ,N. We 

use a grid with N = 500 points in all our simulations. The upper and lower bound of this 
grid are given by cri = —1 + e and (TAr = l —e respectively, where e denotes a UV cutoff. The 
discretized version of our trial solution on this grid is written as 

X^{cji) = X>^ = {Vi,Z^,Xi). (B.l) 


The explicit form of the trial solutions used in the computations of 2-point functions and 
HEE are given in appendix B.l and appendix B.2, respectively. In all our simulations we 
evaluate the cutoff e such that the z-position of the corresponding cutoff surface is fixed at 
Zuv = 0.05, i.e. at rjjv = 20. The boundary conditions are imposed at this cutoff surface 


Vi = VN = t (B.2) 

Zi = Zn = Zuv (B.3) 

Xi = -ll2 Xn = 1/2 (B.4) 

where t denotes the time and I the separation on the cutoff surface. 

The finite difference representation of the geodesic equation (3.2) is given by 

E} = Vi+i - Vi- h{pv)i (B.5a) 

Ef = Zj+i - Zi - h{px)i (B.5b) 

Ef = W+i - Xi- h{px)i (B.5c) 

Ef = {pv)i+i - {pv)i + hJi{pv)i + h{(r^vv)i{pv)‘i + xx)i{Px)j) (B.5d) 

= {Pz)i+i - {pz)i + hJi{pz)i + h{{f^vv)i{pv)‘i 

+2{t^vz)i{Pv)i{pz)i + zz)i{Pz)‘i + xx)i{Px)‘i) (B.5e) 


^i = {Px)i+i — {px)i + hJi{px)i + 2h{{T^xv)i{px)i{pv)i + xz)i{px)i{pz)i) (B.5f) 


where E^ is the residual at point i in equation k] further {pv)i = V denotes the first 
derivatives of the geodesic coordinates; quantitities with bar are averaged like {pv)i = 
. (p . y ) ? .-. i+ . (gy Christoffel symbols xz)i are evaluated from the averaged metric func¬ 
tions; the explicit form of the Jacobian Ji for the case of 2-point functions and HEE are 
given in appendix B.l and appendix B.2. 

The initial guess will in general not satisfy these FDEs very well, i.e. the residua Ef 
will be rather large. To quantify the deviation of a given trial solution to the true solution 
we use the following measure 


6N 


(B.6) 
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The strategy is to compute increments AX^ for the geodesic coordinates such that 
Xnew = ^old IS an improved approximation to the previous solution This we 

do iteratively until we reach 

h < 10-15 . (B.7) 

Equations for the increments are obtained by demanding the hrst order Taylor expansion 
of the hnite difference equations with respect to small changes in the coordinates to vanish. 
We do this following exactly the procedure described in the section on relaxation methods 
in reference [49]. 

The correction AX^ generated from the first order Taylor expansion is in general only 
an improvement close to the true solution. We account for this by introducing a weight a 
that modifies the correction in each relaxation step 

= + (B.8) 

We choose the weight a such that the full correction is used only close to the true solution. 


0.05 if h > 0.05 

1 else 


(B.9) 


In the time evolution we use as ansatz geodesic at time ti = U-i +6t the geodesic form 
the previous time step U-i. With a step size of 6t « 0.05 usually less than 5 relaxation 
steps turned out to be sufficient to reach our accuracy goal of h < 10 “ 


B.l Initialization with Poincare patch AdS geodesics 

The line element in 3-dimensional Poincare patch AdS space is given by 

^ (y—— 2 dz dv + . (B.IO) 

We parametrize the geodesics by an affine parameter r and denote derivatives with respect 
to it by dot. The geodesic equations of motion allow hrst integrals 

X = Lz^ V = Ez'^ — z z = Az\/\ — (L^ — E'^)z'^ (B.ll) 

where L and E are constants of motion and the third equation comes from the spacelike 
condition ds^ = 1. The last equation shows that the geodesics have two branches. It is 
convenient to reparametrize both branches of the geodesics in terms of the holographic 
coordinate z, i.e., to hnd expressions x±{z) and v±{z). We want to initialize the relaxation 
code with geodesics that have a symmetric advanced time with respect to the exchange of 
the branches, V-^-{z) = V-{z). This is achieved by setting E = 0. Then integrating the 
above system yields 

x±{z) = ±\/ L~^ — z"^ v{z) = vq — z (B.12) 

where vq is the boundary time and we have set to zero the additive integration constant in 
x±. In the solution for the advanced time coordinate one clearly sees the ‘bending back in 
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time’-effect as one goes deeper into the bulk. By this we mean simply that v{z) becomes 
more negative as the holographic coordinate z increases. 

With the choices above the constant of motion L is related to the separation I of the 
two endpoints on the boundary 

/ = |x+= |-. (B.13) 

For our numerical algorithm it turns out to be useful to choose a non-afhne parametrization 
that lies in a fixed interval and covers both branches at the same time, namely 

a:((T) = ^(cr\/2 - 0-2) v{a) = vo-z{a) (B.14) 

where a G [—1,1]. We realize a UV-cutoff at a given value Zjjv by truncating the non-afhne 
parater a G [— l-|-e, 1 — e] with e given by 


e = 1 - 



‘^Zjjv 

I 


(B.15) 


To get the Jacobian we need in the non-afhne geodesic equation we hrst integrate the third 
equation in (B.12) to express the affine parameter r in terms of z 


t{z) = ± 


dz 



=Fartanh^y^ 




(B.16) 


Next we use the hrst equation in (B.14) and express the affine parameter r in terms of the 
non-afhne parameter a 

r(cr) = =Fartanh^(T\/2 — . (B.17) 

The Jacobian is then given by 


d^r /dr 5(7 — 3cr^ 
da'^ I da 2 — 3(7^ -|- a'^ 


(B.18) 


B.2 Initialization with conformal AdS geodesics 

We proceed here similarly to appendix B.l. Starting with either of the line elements (2.19) 
in the isotropic static limit {B = 0, T, = r, A = r^) and introducing again z = 1/r yields 


ds^ = ^ (—dr^ — 2d2;d/;-|-dx^) (B.19) 

where we called the third coordinate x, regardless of the specihc case. Integrating the 
geodesic equations, using ds^ = 1 and setting E = 0 yields 


X = Lz^ V = —z z = ±z^^/l — L '^2 


With help of the identity 


dzz^ z^+^ 


Vl - L^z^ 1 + 


* Z7 n l+n 1 I 1+n. r2..61 


(B.20) 


(B.21) 


- 21 















Figure IX. Comparison of the pure AdS geodesic (blue) with the conformal geodesic (red). For 
the same separation the conformal geodesic extends about twice as far into the bulk. 


the third equation in (B.20) can be used to express the affine parameter in terms of the 
holographic coordinate 


T = ± 


dz 

z^\J\ — 



[ 1 , 


1 2 . 
3’ 3’ 



(B.22) 


In order to obtain the solution for the geodesic it turns out to be convenient to solve for x 
in terms of 2 


= ± 


dzLz^ 
\/l — L‘^z^ 


= ^2^ 


Lz^ 

~T~ 


2F1 [ 


1 2 5. r2 6 
2’ 3’ 3’ ^ . 


(B.23) 


The subscript ± corresponds to the two branches of the geodesic and 1/2 is an integra¬ 
tion constant corresponding to the boundary value of the geodesic since the hypergeometric 
function is zero at the boundary z = 0. The endpoint of the geodesic is given hy L'^z^ = 1, 
meaning that the two branches do not join in general. In order to get a smoothly joining 
geodesic we need to adjust L and express it in terms of the boundary separation I 

7rir[5/3f 
8FT [7/6]^ 

Like for the Poincare patch AdS geodesics we choose a non-affine parametrization 

z{a) = Zmax(l - cr^) (B.25) 

x{a) = sgn(c7)(^ - ^ -f 2 F 1 [I, |, ^;L‘^z{af] ) (B.26) 

v{a) = vq — z{a) (B.27) 

where a G [—1,1] and .^max = /^(f) denotes the z-position at which the two 

branches join. We realize a UV-cutoff at a given value Zjjv by truncating the non-affine 


(B.24) 
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parater a G [—1 + e, 1 —e] with e given by 


e = (B.28) 

V ^max 

The affine parameter r in terms of the non-affine parameter a reads 

= ■ (B.29) 

The Jacobian is then given by 

dV /dr _ -51o- + 145^3 - 205o-5 + 159a^ - 65a^ + 
dcr^/ d(T (2 — — (t 2)(3 — Scr^ + (T^)(l — + cr^) 

As can be seen from Fig. IX for the same boundary separation the conformal geodesic 

extends much further into the bulk compared to the Poincare patch AdS case. As a conse¬ 
quence, the boundary separation we can choose for the HEE is much smaller than for the 

2-point functions for any given value of the boundary time. 
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